function [out] = base(A,varargin)
% Finds a base for matrix A. Similar to rref, but faster.
% The 2nd input is the cols to give preference to not dropping. Leave
% empty if no order is important

% The 3rd input is the precision to determine collinearity. Defualt is
% 1e-16 (machine precision).
%#ok<*AGROW> 

[n,m]=size(A); if any(~isfinite(A(:))), error('Matrix contains NaN or Inf'); end
% Define the rcond as Norm * eps * n , where the n is the row being processed
if nargin==3, rcond = sum(abs(diag(A))).*varargin{2};
else rcond = max(abs(diag(A))).*eps; 
end;

% v contains the rows to prioritize on; v2 contains the rest
v=1:n; v2=[]; 
if nargin>1 && ~isempty(varargin{1}), v=varargin{1}; v2=1:n; v2(v)=[]; end;

o=1:n; base=[]; % o is the order of the rows/cols;   
while ~isempty(v), % cycle over priority pivots first, then on non-priority pivots
  % Find col to pivot on (labeled k)
  t=abs(diag(A(v,v))'); tm=max(t); r=length(t); j=max((t==tm).*(1:r)); k=v(j); 
  % Remove vector from list, and restart list if necessary
  v(j)=[]; if isempty(v) && ~isempty(v2), v=v2; v2=[]; end;
  % Columns not yet pivoted
  i=length(base); k2=o(i+1:end); k2(k2==k)=[]; 
  
  if tm>rcond*(i+1), 
   % Grow base and pivot to do only lower triangle reduction
   base=[base k]; o=[base k2];   % Valid base so far, and pivoting
   % Execute row reduction  
   c=[k k2 n+1:m]; r=k2;         % Remaining cols (incl k) and remaining rows (excl k) 
   A(k,c)=A(k,c)./A(k,k); 
   A(r,c)=A(r,c) - A(r,k)*A(k,c); 
  else % Send vector to end; base doesn't change
    o=[base k2 k]; 
  end;
end;

% Return the list of indices of a valid base:
out=sort(base); 

end
